//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#include "ephi.hpp"

static prec_t cali_pot_min, cali_pot_max, cali_pot_add;
static prec_t cali_e_min, cali_e_max, cali_e_add;

static void render (PField& pf, bool first, const std::string& prefix)
{
  std::string str = pf.getChargeConcentrationScript ();
  FILE *fp = fopen ((prefix + "_charge.pov").c_str(), "w");
  fprintf (fp, "camera { location <3, 3, 2>*0.18 look_at <0, 0, 0> sky <0, 0, 1> right x*image_width/image_height }\n");
  fprintf (fp, "#include \"scene.inc\"\n");
  fprintf (fp, "%s", str.c_str());
  fclose (fp);

  str = pf.getPovrayScript ();
  fp = fopen ((prefix + ".pov").c_str(), "w");
  fprintf (fp, "camera { location <3, 3, 2>*0.18 look_at <0, 0, 0> sky <0, 0, 1> right x*image_width/image_height }\n");
  fprintf (fp, "#include \"scene.inc\"\n");
  fprintf (fp, "%s", str.c_str());
  fclose (fp);

  Screen screen(800, 800, 1.0, 1.0);
  screen.set_axis (vect3d(0, 0, 0), vect3d(0, 0, 1).normal(), vect3d(1, 0, 0));
  Statics statics;
  ElectroDynamics dyn(statics);

  PMap pmap(vect3d(0.6, 0.6, 0.6), 0.01, 0.01);
  if ( pf.getMethod() == PField::POINTSYSTEM && 0 )
    {
      pf.fill (pmap);
      pmap.recalc();
      printf ("fill & recalc done!\n");
      statics.addStaticElement (new StaticPMap (pmap));
    }
  else
    {
      statics.addStaticElement (new StaticPotentialField (pf));
    }

  Scene scene(dyn);
  scene.calc(screen, Scene::CALC_POTENTIAL);
  scene.set_coloring (20, false);
  if ( first )
    {
      cali_pot_min = scene.get_fmin();
      cali_pot_max = scene.get_fmax();
      cali_pot_add = scene.get_fadd(); // potential can be negative
    }
  else
    {
      scene.calibrate_field (cali_pot_min, cali_pot_max, cali_pot_add);
    }
  scene.render_map (screen);
  screen.write (prefix + "_pot.bmp");
  scene.calc(screen, Scene::CALC_EFIELD);
  if ( first )
    {
      cali_e_min = scene.get_fmin();
      cali_e_max = scene.get_fmax();
      cali_e_add = scene.get_fadd(); // potential can be negative
    }
  else
    {
      scene.calibrate_field (cali_e_min, cali_e_max, cali_e_add);
    }
  scene.set_coloring (0, false);
  scene.render_map (screen);
  screen.write (prefix + "_efield.bmp");
  scene.set_coloring (20, true);
  scene.render_map (screen);
  screen.write (prefix + "_efield_log_stepped.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  screen.write (prefix + "_efield_lic_logcolored.bmp");
}

int main (int argc, char *argv[])
{
  Ephi::setDebugLevel (Ephi::DEBUG_WARN);
  const char * mnames[] = {"pointsystem", "linesystem", "chargeflow"};
  PField::emethod_t methods[] = {PField::POINTSYSTEM, PField::LINESYSTEM, PField::CHARGEFLOW, (PField::emethod_t)0};
  prec_t potential = 1000.0;
  std::string prefix = "test16_";
  printf ("OK: calculating a torus (should be 31.107pF)\n");
  if (0)
  for ( size_t i = 0; ; i++ )
    {
      if ( !methods[i] )
          break;
      PField pf(methods[i], 0.03);
      pf.addTorus (vect3d(0, 0, 0), vect3d(0, 0, 1), 0.3, 0.075, potential);

      size_t pc, lc;
      pf.getStatistics (pc, lc);
      printf ("OK: created model with %u points, %u links\n", (unsigned int)pc, (unsigned int)lc);

      if ( !pf.equalize() )
        {
          fprintf (stderr, "ERROR: unable to equalize with method %s\n", mnames[i]);
          continue;
        }

      prec_t added_charge = pf.getAddedCharge();
      printf ("OK: method %s capacitance is %.3fpF\n", mnames[i], prec2double(added_charge / potential)*1e12);

      render (pf, i == 0, std::string ("test16_torus_") + mnames[i]);
    }
  printf ("OK: calculating a sphere (should be 27.81625pF)\n");
  for ( size_t i = 0; methods[i]; i++ )
    {
      //PField pf(methods[i], 0.04);
      PField pf(methods[i], 0.02);
      pf.addSphere (vect3d(0, 0, 0), 0.25, potential);

#if 0
      std::vector<vect3d> planes;
      std::vector<transf3d> transforms;
      planes.push_back (vect3d(0, 1, 0));
      planes.push_back (vect3d(1, 0, -1).normal());
      planes.push_back (vect3d(0, -1, 1).normal());
      transforms.reserve (48);
      const transf3d_pair *tr = get_cube_polywell_transforms48 ();
      for ( size_t j = 0; j < 48; j++ )
          transforms.push_back (tr[j].f);
      pf.setTransforms (planes, transforms);

      vect3d n1 = (planes[1] % planes[0]).normal();
      vect3d n2 = (planes[0] % planes[2]).normal();
      vect3d n3 = (planes[2] % planes[1]).normal();

      pf.addSphereSegment (0.25, n1, n2, n3, 1000, false);
#endif

      size_t pc, lc;
      pf.getStatistics (pc, lc);
      printf ("OK: created model with %u points, %u links\n", (unsigned int)pc, (unsigned int)lc);

      if ( !pf.equalize() )
        {
          fprintf (stderr, "ERROR: unable to equalize with method %s\n", mnames[i]);
          continue;
        }

      prec_t added_charge = pf.getAddedCharge();
      printf ("OK: method %s capacitance is %.5fpF\n", mnames[i], prec2double(added_charge / potential)*1e12);

      render (pf, i == 0, std::string ("test16_sphere_") + mnames[i]);
      exit(0);
    }
  return 0;
}

